#ifndef PLANT_SPDE2D_H
#define PLANT_SPDE2D_H

#include "../../include/Signal.h"
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

const double SQRT_TWO = 1.414213562373095;

class plant_SPDE2D:public block
{
  int    mode;
  int    LatticeSize;
  double *Km2[2];

  // States
  double *z1;
  double *z2;

  double T;
  double W;
  double c;
  double c2;
  double sigma;
  double **lambda;

  double dt_con; // dt_con is the dt used in numerical integration. 
  // block::dt is the time interval for output.
  double sqrt_dt_con;
  double dt_conservative_factor;

  // Random number generator supported by GSL
  const gsl_rng_type * noise_type;
  gsl_rng * noise;

  double eq_c(double W);
  double eq_c2(double W);
  double eq_sigma2(double W);

  inline int index1D(int m){return m;};

  inline double phi1m(int m,int ix){
    return sqrt(2/Xmax)*sin(2*m*pi*ix/Xmax);
  };
  //double phi2mn(int m,int n,int ix,int iy);
  inline double phi2m(int m,int ix){
    double coeff;
    if((m == 0 ){
      coeff = sqrt(1/Xmax);
    }
    else{
      coeff = sqrt(2/Xmax)*cos(2*m*pi*ix/Xmax);
    }
    return coeff;
  };


public:
  plant_SPDE2D(int M,int L,double idt);
  ~plant_SPDE2D();
  virtual void update(double sysTime);
  virtual void reset();
  void update(double sysTime,double T,double W);
  void output(string filename);
  double getM2();
  double getR2();
};

void h2mode_2D(double *h,int L,double *z,int M);
void mode2h_2D(double *z,int M,double *h,int L);
double randn();
#endif